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Abstract 

We compute by Monte Carlo methods the individual distributions of the kth smallest Dirac 
operator eigenvalues in QCD, and compare them with recent analytical predictions. We do this for 
both massless and massive quarks in an SU(3) gauge theory with staggered fermions. Very precise 
agreement is found in all cases. As a simple by-product we also extract the microscopic spectral 
density of the Dirac operator in SU(3) gauge theory with dynamical massive fermions for Nf = 1 
and 2, and obtain high-accuracy agreement with analytical expressions. 
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In the chiral limit of QCD it is possible to derive exact analytical expressions for eigenvalue distribu- 
tions and eigenvalue correlations of the Dirac operator. This situation arises if one restricts oneself 
to the phase of spontaneously broken chiral symmetry (with infinite- volume zero-mass condensate S), 
and considers the theory in a finite space-time volume V. Taking the chiral limit by sending quark 
masses rrii to zero in such a way that the combination /ij = uiiTiV is kept fixed as ^ ^ oo, this leads to 
an interesting finite- volume scaling regime of QCD. Denoting generically all pseudo-Goldstone masses 
by m^, exact scaling is achieved if m,r V^~^/^ The only assumption is that chiral symme- 
try is spontaneously broken from one group G to a smaller group H. The coset G/H specifies the 
universality class [Q]. 

The most detailed analytical predictions were originally obtained from an intriguing exact relation to 
universal Random Matrix Theory results ||2|, Later, the same results have been derived directly 
from the effective partition functions Q based on the so-called supersymmetric formulation |p. Very 
recently, also the replica method Q has been shown to yield identical results ||^]. Beyond any doubt 
now, these analytical predictions represent exact statements about the Dirac operator spectrum of 
QCD in the above finite- volume scaling region. This is an exciting situation from the point of view of 
lattice gauge theory, which is almost tailored to study such a finite-volume regime. It is one of the few 
instances where there are exact and non-perturbative analytical predictions that immediately can be 
compared with lattice Monte Carlo results. There have by now been several Monte Carlo studies of 
these analytical predictions for Dirac operator eigenvalues in various four-dimensional gauge theories 
m, ^ 10, 11, |l^, 13 1 . One peculiarity of these predictions is that they are most easily expressed in 



gauge field sectors of fixed topological charge i^. For fermions sensitive to topology this gives an 
interesting new possibility, as predictions differ markedly for different topological sectors [0]. But 
staggered fermions seem oblivious to gauge field topology at the gauge couplings that are presently 
realistic |15, predictions can then only be made with respect to the topologically trivial sector of 



1^ = 0. 

In this letter we shall extend the Monte Carlo analyses of Dirac operator eigenvalues in two directions. 
First, we shall confront some very recent analytical results concerning the probability distribution of 
the kth smallest Dirac operator eigenvalue with Monte Carlo data. This is important, because 
there is much more detailed information in these individual eigenvalue distributions than in, for ex- 
ample, the spectral density, which is the sum of all these individual distributions (see below). Second, 



we shall make the extension of the quenched results of Ref. [10| to the case of dynamical fermions of 
non-zero mass. Until now, Monte Carlo analyses of microscopic Dirac operator spectra with dynamical 
fermions have been limited to a different universality classQ than the one relevant for QCD Q (see also 
Ref. |l^]). We study gauge group SU(3), which, even with staggered fermions, allows us to compare 
with the universality class relevant for continuum QCD (the class of Dyson index f3 = 2). 



Details of how we compute the lowest-lying Dirac operator eigenvalues have been given in Ref. [|1| 
and we shall not repeat them here. Instead, let us briefly recall the analytical predictions with which 
we should compare our Monte Carlo data. We start with the recently derived expression for the 
distribution of the feth Dirac operator eigenvalue in the finite-volume scaling region mentioned above. 
The generator for these individual eigenvalue distributions is the joint probability distribution for the k 
smallest eigenvalues. It turns out that this joint probability distribution can be written very compactly 
in terms of finite- volume partition functions with additional fermion species. For the universality class 

^The sum over topology can be performed [|l^, but closed analytical forms are not known except in very simple cases. 
For instance, in theories with massless (dynamical) fermions only the sector with i/ = contributes. 
^An SU(2) gauge theory labeled by the so-called "Dyson index" (3 = 4 for staggered fermions H. 
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relevant for QCD, the explicit expression is |17] 
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Here the additional 2(fc — 1) fermion species of masses yCl ~ Q ^ire doubly degenerate, as indicated. 
The partition function in the numerator is evaluated in a sector of fixed topological charge v = 2 
(independent of the actual topological charge in question). But there is information about v in the 
fact that V additional fermion species (all degenerate in mass Cfc) enter in the partition function. 
Actually, since we shall be dealing with staggered fermions in the present paper, only the sector with 
J/ = is relevant. In that case the above formula simplifies: 

fc-i fc-i 



2o({/x}) 



(2) 



The fact that in particular the smallest (non-zero) Dirac operator eigenvalue has a distribution deter- 
mined entirely in terms of the effective partition function was noted earlier , but the present most 
general expression is actually more compact. The proportionality factor C depends on the chosen 
normalization of the partition functions. We choose a normalization in which [20, Q| 

^,({m}) = det^({/x})/A(^2) ^ (3) 

where the determinant is taken over the N j x A'^^ matrix 

^({/^D = f4''lu+j-M , (4) 

and 

Nf 

a(m') = n(/^'-^i)- (5) 

i>j 

With this convention the normalization factor is C = 1/2 for all values of k, Nf and ly. 

For the purpose of comparing with lattice gauge theory data, a more convenient quantity to focus on 



is the probability distribution of the kth. smallest Dirac operator eigenvalue |17|. One gets it from the 
joint probability distribution by integrating out the previous k — 1 smaller eigenvalues: 
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The individual eigenvalue distributions by definition sum up to the microscopic spectral density 

oo 

ps{Q■,{^A) = Y.Pk{C;{p})- (7) 

k=l 



2 



0.5 




Figure 1: Simulation data of the 10 lowest quenched Dirac operator eigenvalues compared with the 
analytic predictions for the individual eigenvalues. The data were obtained on an 8^ lattice at (3 = 5.1. 



But a much simpler way to get the microscopic spectral density is by direct evaluation. As we will in 
addition be interested in the case of dynamical fermions with non-zero quark masses /ij, we need the 



microscopic spectral density for the general massive case of this universality class 22, 23 1. Also 



this quantity can most simply be expressed in terms of the finite- volume partition functions |Q] : 

Nf 

psic-Af^}) = (-incin(c'+M|)^4^#%^ , (8) 

using the same normalization convention as above. This spectral density is "double-microscopic" in 
the sense that both eigenvalues A and masses rrij have been rescaled according to C = XTiV and 
jij = rrijTiV . As the rescaled masses p,j are sent to infinity, one obtains a series of decoupling relations 
where the corresponding quark flavors become effectively quenched [^]. Because we precisely want 
to see the effects of having dynamical fermions, it is thus important to keep the rescaled masses //j on 
roughly the same scale as the eigenvalues, or smaller. 



We begin by comparing quenched lattice gauge theory data with the predictions for the individual 
eigenvalue distributions. We have here large statistics for the 10 smallest Dirac operator eigenvalues 
on an 8^ lattice at coupling j3 = 5.1. These parameters and all others used in this paper are summarized 
m Table ^ In Fig. | we show a comparison of these lattice Monte Carlo results with the formula (y) 
for = 1 . . . 10. S was determined from a fit of the distribution of the lowest eigenvalue. The formula 
obviously gets quite tedious to evaluate by numerical integration for large values of k. The analytic 
curves for the distributions of the last couple of eigenvalues have actually not been obtained directly 
from (^), but rather from the eigenvalue distributions of computer generated random matrices (which 
obviously amounts to the same). The agreement with our lattice Monte Carlo data is nothing less 
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V 


N/ 


13 




m 


#configs 


44 


1 


4.6 


1.24 


0.02 
0.01 
0.004 


40,600 
32,650 
25,800 


44 


2 


4.2 


1.28 


0.01 


26,040 


64 


1 


4.7 


1.17 


0.003 


3,370 


84 





5.1 


1.15 




17,454 



Table 1: Simulation parameters. Note that the infinite- volume zero-mass chiral condensate S is one 
single free parameter, depending only on j3 and Nj. For the case of y = 4"^ and Nf = 1 there are thus 
no free parameters at all for two of three simulations with different quark masses. 

than spectacular for all the eigenvalues we have available. Eventually, on any given finite lattice, 
1/y-corrections will take over, and this very precise agreement with the theoretical predictions will 
no longer hold. But in the limit V ^ oo the agreement should persist for an infinity of smallest 
Dirac operator eigenvalues A^, rescaled so that Qn = AnSy remains finite in the infinite-volume limit. 
Agreement with these lattice data and the microscopic spectral density has already been presented 



in Ref. |15|, but after just a few oscillations the spectral density quickly becomes an almost constant 
function, with very little information. For the first time we see here that there is detailed agreement 
with the analytical predictions for a whole sequence of individual eigenvalues underneath. 

We next turn to simulations with dynamical fermions of finite mass. As mentioned above, there has 
previously been only one Monte Carlo study in this direction for a different gauge group (and 
hence different universality class) SU(2). We shall here provide first results for a study with gauge 
group SU(3), and hence the universality class of continuum QCD. We shall not only compare with 
the microscopic spectral density, but also focus on individual eigenvalue distributions, just as we did 
above in the quenched case. 

For our dynamical simulations we employed the Hybrid Monte Carlo algorithm with Nf = 1 or 2 
species of staggered fermions. In the continuum limit, this would correspond to nj = ANf continuum 
fermion flavors. However, in order to have a large physical volume with few lattice sites — our 
computer resources limit us to the use of systems with 44 up to 64 lattice sites and the finite size 
scaling limit we want to consider requires a sufficiently large physical volume, V ^ ^qcd — 
forced to work at a large lattice spacing and hence a large gauge coupling. In such a situation the 
flavor symmetry of the staggered fermions is badly broken, especially in the limit of small masses 
of interest here. Only the Nj true pseudo-Goldstone pions from the spontaneous breaking of the 
U{Nf) X U{Nf) symmetry to the diagonal U{Nf) symmetry are light. The flavor symmetry breaking 
makes the other — 1 — "would-be" pseudo-Goldstone pions much heavier. As far as the infrared 
limit is concerned, only Nf flavors contribute. This will be verified by our results presented below. 

In Fig. I we compare the spectral density of the simulation data for different quark masses. S was 
determined, for each of the three ensembles, from a fit of the distribution of the lowest eigenvalue. 
The values agreed within errors, and were averaged to obtain the result in Table ||. This agreement 
is of course as it should be: S is the infinite-volume chiral condensate in the chiral limit m — s- (and 
in particular mass- independent). For each mass we plot the analytical result (^ for that particular 
mass value (rescaled according to fi = mT,V), and include the curve for the limiting cases of the 
quenched result corresponding to infinite mass quarks, as well as for one massless quark. We see 
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Figure 2: The microscopic spectral density for different quark masses on a 4^ lattice with Nf = 1 
at P = 4.6. The three simulation masses m = 0.02,0.01 and 0.004 correspond to /i = 6.35,3.17 and 
1.27, respectively. Also the curve for one massless flavor and the Nj=0 (quenched) curve, which is 
equivalent to the limit fi — > oo, is plotted. There is excellent agreement with the analytical curves 
corresponding to just the right ^-value. 




Figure 3: The smallest eigenvalue distributions of the same data as in Fig. |2[ Agreement holds also 
at individual distributions. The sum of the individual eigenvalues, the microscopic spectral density, is 
shown by the upper curve on all three plots. 



that the data beautifully interpolate between the two limits and follow the analytical curve almost 
perfectly. Looking at the individual eigenvalue distributions in Fig. ^ we get a better view on how far 
the agreement goes. It is clear that only the very first eigenvalues match perfectly and that there is 
a very tiny discrepancy for the last (of the 6) distributions due to the limited volume, but still the 
agreement is impressive. We did a single run on a larger volume (6*^), and indeed the agreement is 
better, see Fig. Q B. Finally, data from a run corresponding to two flavors of quarks are presented in 
Fig. ^ A, with the same degree of agreement. 

To conclude, we have found very precise agreement with the analytical predictions for the distributions 
of Dirac operator eigenvalues in the microscopic finite-volume scaling region. This is seen in both 
quenched and unquenched simulations, and for the first time there are now high-statistics comparisons 
between the analytical predictions for the Dirac operator spectrum of massive quarks and lattice Monte 
Carlo data for gauge group SU(3). Our Monte Carlo calculations have been restricted to staggered 
fermions, which are insensitive to gauge field topology at the couplings considered here. We have 
thus compared with analytical predictions for the topological sector of charge z/ = 0. Perhaps the 
biggest surprise is that the previously seen agreement with analytical predictions for the summed- 
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Figure 4: The smallest eigenvalue distributions for (A) the Nf = 2 theory with degenerate masses on 
a 4"^ lattice, and (B) the Nj = 1 theory on a 6^ lattice. Gauge couplings as shown in Table 0. Again 
the spectral density corresponding to these Nf and fj, parameters are shown as the upper curves. 

up spectral density holds at a much more detailed level underneath. The individual smallest Dirac 
operator eigenvalues are really completely locked up by gauge field topology, and, away from the exact 
zero modes, by the coset G/H oi spontaneous chiral symmetry breaking. 
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